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WSTĘP 


W ostatnich kilkunastu latach opracowano kilka numerycznych metod 
symulujących procesy elektrochemiczne. Są to przede wszystkim metody różnic 
skończonych (jawna i uwikłana) [1, 2] oraz opracowana w latach 70. metoda kol- 
lokacji ortogonalnej z użyciem rozwinięć symulowanych funkcji w szereg wielo- 
mianów ortogonalnych [3]. Najbardziej prostą z nich jest metoda jawnych różnic 
skończonych opracowana przez Feldberga [1]. Prater i Bard zaadaptowali tę metodę 
do symulacji procesów zachodzących na wirującej elektrodzie dysk-pierścień. 

Celem niniejszej pracy jest szczegółowe opisanie metody Pratera i Barda 
oraz prezentacja zasadniczych części napisanego w języku Fortran programu do 
symulacji procesu, w którym substancja generowana na dysku podczas ruchu w 
kierunku pierścienia ulega homogenicznej reakcji w warstwie dyfuzyjnej zgodnie z 
kinetyką pierwszego rzędu (proces EC I). Do tej pory nie publikowano nigdzie pro- 
gramu do symulacji procesów zachodzących na wirującej elektrodzie dysk-pierścień. 
Metoda ta ma kilka podstawowych zalet. Jest bardzo prosta, łatwo daje się przy- 
stosować do badania bardziej złożonych procesów np. procesów EC II rzędu [5], 
procesów ECE [6] itp. Ponadto, jak każda z metod symulacyjnych, może być sto- 
sowana bez ograniczeń związanych z wymiarami elektrody. Jak wiadomo równania 
analityczne opisujące rozkład stężenia pod elektrodą dysk-pierścień mają charak- 
ter przybliżony; stosuje się je przede wszystkim w przypadku elektrod z cienkim 
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pierścieniem i cienką przerwą izolacyjną (thin ring-thin gap) [7]. Natomiast w prak- 
tyce trudno jest skonstruować elektrodę spełniającą powyższy warunek, dlatego 
lepiej jest prowadzić badania wykorzystując metodę symulacji komputerowych. 


DYSKRETYZACJA PRZESTRZENI POD ELEKTRODĄ DYSK-PIERSCIEN|4] 


Większość problemów elektrochemii została sformułowana za pomocą pojęcia 
kontinuum i daje się opisać za pomocą układów równań różniczkowych. Jednakże w 
przypadku rozwiązywania problemów elektrochemicznych przy wykorzystaniu sy- 
mulacji cyfrowej sachodsi potrzeba stosowania modeli dyskretnych i skończonych. 
Wymaga to zastąpienia kontinuum siatką przestrsenno-czasową. 

Obszar pod elektrodą został podzielony na pewną liczbę cienkich, równoleg- 
łych warstw (o grubości Ar cm). Elektroda dysk-pierścień jest umieszczona w 
środku pierwszej warstwy [8]. Z kolei każda warstwa podzielona jest w kierunku 
radialnym na jeden duży cylinder o promieniu równym promieniowi dysku i środku 
w osi rotacji (cylinder ten to obszar leżący pod dyskiem; takie rozwiąsanie podziału 
przestrzeni w kierunku radialnym wynika z założenia, że powierzchnia dysku jest 
jednakowo dostępna dla substancji) oras szereg koncentrycznych pierścieni o sze- 
rokości Ar. Przedział Ar musi spełniać warunek: 


112 (IRI = 0.5)Ar 


gdzie rı jest promieniem dysku w cm. IRI jest promieniem dysku wyrażonym w 
jednostkach Ar. Powstaje w ten sposób siatka objętościowych komórek. Każdej z 
nich przyporządkowuje się parę liczb (J, K). J wzrasta w kierunku normalnym 
do elektrody, a K w kierunku radialnym. Duża cylindryczna komórka w pierwszej 
warstwie (J = 1, K = 1) odpowiada elektrodzie dyskowej, natomiast koncentryczne 
pierścienie w warstwie pierwszej (J = 1) odpowiadają, kolejno, przerwie izolacyjnej 
i elektrodzie pierścieniowej. 


OGÓLNY SCHEMAT SYMULACJI 


W chwili t = 0 stężenia w poszczególnych komórkach wynikają z przyjętych 
dla danego procesu warunków początkowych. Symulacja rozpoczyna się smianą 
warunków początkowych na powierzchni elektrody dyskowej (komórka (1, 1)) sgo- 
dnie z przyjętymi warunkami brsegowymi. Zmiana ta, odpowiadająca prsedsialowi 
At, wymusza powstanie gradientu stężenia pomiędzy powierzchnią elektrody a ros- 
tworem. Zmiany stężenia w komórkach swiąsane z dyfuzją w kierunku normalnym 
(dyfuzja radialna jako mała jest pomijana [9]) są oblicsane prsy użyciu równań 
różnicowych będących dyskretnym odpowiednikiem równań różniczkowych Ficka. 
W celu obliczenia smian stężenia spowodowanych efektami hydrodynamicsnymi 
należy określić odległość jaką przebył rostwór w czasie At w kierunku radialnym i 
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normalnym. Następnie modyfikuje się odpowiednio stężenia w komórkach, tsn. jeśli 
roztwór znajdujący się w komórce (J, K) w chwili t+ At snajdował się w chwili t w 
komórce (J + AJ, K- AK), to należy zastąpić stężenie w komórce (J, K) stężeniem 
s komórki (J + AJ, K — AK) [4]. Kinetykę reakcji homogenicznej uwzględnia się 
obliczając dla każdej komórki ilość substancji, jaka w niej przereagowała w czasie 
At wykorzystując prawo szybkości reakcji. Po obliczeniu efektów dyfuzyjnych, kon- 
wekcyjnych i kinetycznych otrzymujemy rozkład stężenia w komórkach po czasie 
At (jedna iteracja). Należy następnie smodyfikować stężenia na powierzchni całej 
elektrody zgodnie s warunkami brzegowymi i przejść do kolejnej iteracji dla 2At. 
Iteracje są realizowane do osiągnięcia stanu stacjonarnego. Symulacja będzie tym 
dokladniejssa im mniejsze będą przedziały At, Ar i Az. 


WARUNKI POCZĄTKOWE I BRZEGOWE 


Rospatrywany proces można przedstawić prsy pomocy reakcji: 


Atne—B (1) 
B— X (2) 
B ne A (3) 


W chwili t = 0 potencjał na elektrodzie dyskowej osiąga potencjał prądu 
granicznego. Potencjał pierścienia jest taki, aby cała ilość substancji B docierająca 
do pierścienia była natychmiast przekształcana w A. Ponieważ wszystkie obliczenia 
w symulacji operują wielkościami beswymiarowymi wygodnie jest używać stężeń 
odniesionych do stężenia substancji A. I tak [4]: 


FA=C,/C3, TB OO 


gdsie FA i FB są tablicami reprezentującymi stężenia substancji A i B w rostworse, 
O jest stężeniem substancji A w chwili t = 0, a C4 i Cg są stężeniami A i B w 
danym momencie procesu. W chwili t = 0 mamy następujące stężenia substancji 
biorących udsiał w procesie: 

FAI1(J, K) = 1.0 


FBI(J, K) = 0.0 
Po smianie potencjału na dysku do potencjału prądu granicznego mamy następujące 
warunki brzegowe: 

FAi(1, 1) = 0.0 

FBi(1,1) = 1.0 
(czyli cała ilość substancji A docierająca do dysku ulega reakcji elektrodowej do 
substancji B) 
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oras 


FAi(1, K,) = 1.0 
FB1(1, K,) = 0.0 


(csyli cała ilość substancji B docierająca do pierścienia ulega reakcji elektrodowej 
do substancji A, przy csym Kp oznacza wartości K odpowiadające elektrodzie 
pierścieniowej). 


DYFUZJA 


W wyniku zmiany potencjału powstaje gradient stężenia pod dyskiem 
pomiędzy pierwszą (J = 1) a drugą (J = 2) warstwą, co powoduje dyfuzję substan- 
cji do elektrody. Wykorzystując drugie prawo Ficka w wersji różnicowej możemy 
napisać [1]: 


FA(J, K) = FAI, K) + DMA(FA1(J + 1, K) — 2FA1(J, K) + FA1(J — 1, K)) 
FB(J, K) = FB1(J, K) + DMB(FBI(J + 1, K) - 2FBI(J, K) + FBi(J — 1, K)) 
gdzie tablice F A i FB reprezentują stężenia w komórkach po czasie At. Bezwymia- 
rowe współczynniki DMA i DMB są równe [1,4]: 

DMA = D,At/(Az)?; DMB = DpgAt/(Az)? 
gdzie DA i Dg są współczynnikami dyfuzji substancji A i B. 


Współczynniki DMA i DMB muszą być mniejsse od 0.5 [1], aby proces 
symulacji był zbieżny. Najczęściej zakłada się, że DM A = DMB (choć założenie to 
nie jest konieczne). 


KONWEKCJA 


Jeśli snana jest prędkość przepływu rostworu (spowodowana obrotem elek- 
trody) dla danej komórki (J, K) w chwili t, łatwo jest obliczyć położenie rostworu 
w chwili t At. W pobliżu elektrody ruch rostworu w kierunku normalnym można 
opisać wyrażeniem [7]: 

v. = 0.51% % 4/222 (4) 
gdzie w jest szybkością obrotów elektrody w rad/s, v jest lepkością kinematyczną 
rostworu w cm?/s, z jest odległością od elektrody, a v, jest prędkością rostworu w 
kierunku z. Ponieważ v, = dz/dt mamy 


dz/dt = —0.5103/2,71/272 (5) 


Po rozwiązaniu równania różniczkowego otrzymujemy 


1/42 — 1/24 = 0.510% % 0 — ti) (6) 


gdzie cą jest odległością rostworu od elektrody w chwili tą a 21 jest odległością w 
chwili ti. Odległość z2 można wyrasić prsy pomocy parametru J, mianowicie [4]: 


z2 = (J — AZ = XJ J Az (7) 


gdzie XJJ osnacsa odległość od elektrody wyrażoną w jednostkach Az w chwili 
tą. Wartość cy można wyrazić wzorem: 


z, MJA (8) 


gdzie XJ jest odległością roztworu od elektrody w chwili ty. Wstawiając wyrażenia 
(7) i (8) do (6) otrzymamy: 


1/XJJ — 1/XJ =0.51w*/3v7124ĄzAt (9) 


lub 
XJ =XJJ/[1 - XJJ(0.510*/2v7W24zAt)] (10) 


Wyrażenie w nawiasie okrągłym oznaczamy V ZERO |4]: 
XJ = XJJ/(1— XJJ VZERO) (11) 


Podobnie postępujemy w prsypadku konwekcji radialnej rostworu. W pobliżu elek- 
trody mamy [7]: 
v, = 0.51w*/3y71/2gy (12) 


gdzie v, jest prędkością rostworu w kierunku radialnym, r jest odległością od osi 
rotacji. Ponieważ v, = dr/dt, więc można zapisać: 


dr / dt O. 5 10% 7 aer (13) 
i po roswiąsaniu równania: 
In(r;/ra) = —0.51v3/2y71/2g(tą — ty) (14) 


Tı i ra oznaczają odległości od osi rotacji rostworu w chwilach, odpowiednio, ti i 
tą. Podstawiając sa z otrzymamy: 


In(r1/ ra) = —0.51w*/3v7/2X JAzAt = -V ZERO XJ (15) 

Jeśli przedstawimy r; i ra prsy pomocy parametru K to otrsymamy [4]: 
rg = |(K + IR1) - 2]Ar = RKAr (16) 
r, = RKKÓr (17) 
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i ostatecznie: 


RKK = RKexp(-VZERO XJ) (18) 


gdzie RK i RKK oznaczają odległości roztworu od osi rotacji w jednostkach Ar 
w czasie, odpowiednio, tą i tı. Para wartości XJ i RKK określa więc położenie 
roztworu w chwili ti, w sytuacji gdy w chwili tą znalazł się on w komórce (J, K). 
Należy więc aktualne stężenie w komórce (J, K) zastąpić stężeniem roztworu o 
położeniu (XJ, RK K), przy czym odległości XJ i RKK musimy przeliczyć na nu- 
mer komórki odpowiadającej temu położeniu. Ilustruje to schemat: 


inn i irr oznaczają numer 
komórki obliczony z wartości XJ i RKK 
Modyfikacja przeprowadzana jest w następujący sposób: 


FAI1(J, K) = FA(INN, IRR) + P (20) 


BIG, K) = FB(INN, IRR) + P (21) 


gdzie tablice FA1 i FBI reprezentują stężenia w komórkach po uwzględnieniu 
konwekcji a P jest „poprawką” zwiąsaną z tym, że wartość stężenia w komórce 
(INN, IRR) odpowiada stężeniu w centrum komórki, a punkt o współrzędnych 
(XJ, RKK) może leżeć na jej obrzeżach (tak jak to pokazano na schemacie). Po- 
prawki liczone są przy założeniu liniowości zmian stężenia pomiędzy sąsiednimi 
komórkami [1]. Ponieważ w przypadku wirującej elektrody mamy do czynienia 
z procesem stacjonarnym (tzn. v, i vz nie zmieniają się w czasie), to obliczenia 
związane z ustaleniem wartości XJ i RKK (a co za tym idzie także wartości INN 
i IRR) dla wszystkich komórek (J, K) wystarczy wykonać raz na początku danej 
symulacji. 

Parametr V ZERO jest bardzo ważny, gdyż wraz z parametrem DMA ustala 
wielkość przedziałów Az i At. Ponieważ Ax = (D4 At/DMA)*/2, to 


V ZERO =0.51v3/2D)/%v-/2DMA7"V/2At*/2 (22) 


Wprowadzając parametr L oraz ty, które są związane zależnością tę = LAt [4], i 
podstawiając tL“! za At w wyrażeniu (22) otrzymamy: 


V ZERO = 0.510% D}? v- 1/2 DM A" 1248/2 L-3? (23) 
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Zakładając, że tę = 9 17 p7" (0.51% (wyrażenie po prawej stronie ma 
wymiar czasu) mamy: f 
V ZERO = DMA 7½L—59 24) 


Ustalając więc arbitralnie jakąś wartość DMA (najczęściej DMA > 0.4, co daje 
szybkosbieżną symulację) możemy poprses manipulowanie parametrem L zwiększać 
lub smniejssać dokładność symulacji. Im większa wartość L tym mniejsse mamy 
prsedsialy At i Az, co wpływa na dokładność obliczeń, zwiękssa jednak liczbę 
wykonywanych iteracji. 


KINETYKA 


W przypadku kinetyki pierwszego rzędu (reakcja (2)) mamy sgodnie s pra- 
wem szybkości reakcji: 
dcp / dt = —keg (25) 


co można zapisać w wersji różnicowej używając wielkości z symulacji: 
AFBi(J, K) = —kAtFBi(J, K) (26) 
lub podstawiając za At [4]: 
AFBA(J, K) = -|kw-1v7-1/5D71/3(0.51)-2/3]FB1(J, K)/L (27) 
Wyrażenie w nawiasach kwadratowych oznaczamy X KT: 
AFBI(J, K) = —XKT FBi(J,K) (28) 


lub 
FB1(J, K) = FB1(J, K) (I- XKT/L) (29) 


ponieważ FB1(J, K) = FB1(J, K) — AFBI(J, K). 


OBLICZENIA WARTOŚCI PRĄDU 


W obliczeniach natężenia prądu na dysku wykorzystuje się beswymiarowy 
parametr ZD [4]: 


ZD =ip/(0.51)1/3nFApcy DZ?w'/2y-1/6 (30) 


gdzie Ap oznacza powierzchnię dysku a tp jest natężeniem prądu na dysku. Para- 
metr ZD można wyrasić prsy pomocy wielkości s symulacji: 


ZD = DMA'/2FA1(2,1)L'/2 (31) 
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Podobnie, natężenie prądu na pierścieniu wyraża się beswymiarowym parametrem 
ZR [4]. Mianowicie: 
ZR = ip/(0.51)?/*nF Ape®, DZPw!/2y71/6 (32) 


lub 
ZR = L?DMB"/? 5 > FB1(2, K)A(K)/Ap = 
K, 


L'2DMBI? FBI, K)|[(K+IR1—1.5)?- (K+1R1—2.5)?]/(TR1—0.5)? (33) 
K, 
gdsie Kp oznacza wartości K dla pierścienia, A(K) jest powierschnią K-tego 


pierścienia a Ap jest powierzchnią dysku. Stosunek ZR do ZD jest efektywnością 
sbierania elektrody dysk-pierścień N. 


MODYFIKACJA WARUNKÓW BRZEGOWYCH 


Po oblicseniu efektów swiąsanych s dyfusja, konwekcją i kinetyką reakcji 
homogenicznej należy smodyfikować warunki brsegowe, co kończy pierwszą iterację. 
Na dysku mamy [4]: 

FA\(1, 1) = 0.0 
FB1(1, 1) = FB1(1, 1) - DMB BI, 1) - FBI, i)) + DMA FA1(2, 1) 
W przerwie izolacyjnej: 
FA1(1, K) = FA1(1, K) + DMA(FA1(2, K) — FA1(1, X)) 


FB1(1, K) = FBi(1, K) + DMB(FB1(2, K) — FBi(1, K)) 
gdzie wartości K odpowiadają przerwie isolacyjnej. 
Na pierścieniu: 
FA1(1, K) = FA1(1, X) + DMB FBi(2, K) - DMA(FAi(1, K) — FA1(2, X) 


FB1(1,K) = 0.0 


gdzie wartości K odpowiadają elektrodzie pierścieniowej. Po smodyfikowaniu wa- 
runków brsegowych następuje wykonanie kolejnej iteracji dla t = 2At, następnie 
dla t = 3At itd. aż do osiągnięcia stanu stacjonarnego. Zastępowanie wartości 
„starych” stężeń (na początku prsedsialu At) stężeniami smodyfikowanymi (na 
końcu prsedsialu At) w obrębie jednej iteracji odbywa się sgodnie se schematem: 


warunki początkowe 
t=0 
TAB1(J,K) 


modyfikacja 
c warunków 
t >0 —— pocratkowych 


dyfuzja 
TAB(J,K) -> TAB1(J,K) J=1,TAB1(1,K) -> TAB1(1,K) 


konwekcja 
TAB(J,K) -> TAB1(J,K) 


kinbtyka 
TAB1(J,K) -> T4310.) öͤͥ⁊uau 
PROGRAM 

Program do symulacji cyfrowej został napisany w języku Fortran. Poniżej 
prezentowane są zasadnicze części programu z krótkim omówieniem. Pominięte so- 
stały te części programu, które związane są z operacjami wejścia-wyjścia, takimi jak 
komunikacja s użytkownikiem programu, tworzenie zbiorów danych, wyświetlanie 
wartości kontrolnych itp. 


1) Ustalenie parametrów początkowych procesu 


dma = 0.46 

dmb = 0.46 

iri = 

ri = 

r2 = 

r3 = 

xkt = 

eps = 

1 = 60 

vsero = (1./eqrt(dma))*(1**(-1.6)) 
delr = ri/(iri-.6) 

ir2 = r2/delr*1 

ir3 = r3/delrei. 

ki = ir2-iri+i 

k2 = ir3-iri+i 

lx = 6*(dma/vzero)**.33333333333+.6 


gdsie rl, r2, r3 oznaczają rzeczywiste wymiary elektrody dysk-pierścień 
(można zamiast tych wymiarów podawać wymiary zredukowane r; = r,/r;, rż = 
12/r1, r3 ri), xkt = kw”! D"!/3u1/3(0.51)72/3 jest parametrem związanym se 
stałą szybkości reakcji k; jeżeli symulowany proces jest opisany reakcjami (1) i (3) 
to xkt = 0.0. Wartość 1 za Praterem i Bardem ustalono na 50, co pozwala skrócić 
czas wykonania programu, przy czym konieczne jest zastosowanie współczynników 
korekcyjnych w równaniach (11) i (18) (patrz punkt 3), oras wartości dma i dmb 
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muszą być równe 0.45. Lx jest ilością warstw w kierunku normalnym do elektrody; 
wartość tego parametru wynika z kombinacji wyrażeń [1]: 


= 1.61(D/v)'/3(v/w)'/2 


70 =0.51(v0*/v)1/2 
gdsie 6 jest grubością warstwy dyfuzyjnej a r" jest odpowiednikiem beswymiarowej 
wielkości vsero wyrażonym w (cm s)". Delr jest wartością Ar, irl, ir2, ir3 są wy- 
miarami elektrody wyrażonymi w jednostkach Ar. K1 i k2 określają ilość komórek 
w kierunku radialnym przypadających, odpowiednio, na dysk s przerwą izolacyjną 
oras na dysk s przerwą i pierścieniem. Parametr irl należy dobrać tak, aby były 


spełnione warunki: 
(ir2-0.5)/(ir1-0.5) = r2/r; oraz (ir3-0.5)/(ir1-0.5) 1/1 


2) Warunki początkowe 


8) Ustalenie wartości zj i rkk związanych z hydrodynamiką procesu 


in(1) = 1 
w(1) = 0. 
do 30 j = 1,1x 
xjj = float(j) - 1. 
xj = xjj/(i. - 1.11 + xjj * vzero) 
in(j) = xj + 1 
w(j) = xj - float(in(j)) + 1. 
do 30 k = 1,k2 
rk = (float(k + iri) - 2.) 
rkk = rk * exp(-1.03 * vzero * xj) 
ir(j,k) = rkk - iri + 2 
30  wi(j,k) = rkk - float(ir(j,k) + iri) + 2 


Wartości xj i rkk odpowiadające położeniu rostworu w chwili t; są prseli- 
czane na numery komórek i sapamiętane w tablicach in(j) i rk(j,k). Ponieważ punkty 
o współrzędnych (xj,rkk) snajdują się najczęściej w obrębie komórki posa jej cen- 
trum, tablice w(j) i w1(j,k) sapamiętują położenie rostworu w obrębie komórki w 
stosunku do jej środka (patrz schemat). Bedsie to potrzebne pray oblicsaniu popra- 
wek w czasie ustalania smian stężeń w komórkach spowodowanych ruchem rostworu 
pod obracającą się elektrodą. 
xjj - odległość od elektrody rostworu s warstwy j w chwili ta, 
rk - odległość od osi rotacji rostworu s komórki (j,k) w chwili ta. 
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We wzorach na xj i rkk użyto czynników korekcyjnych związanych z ustale- 
niem wartości parametru | na 50, które powodują skrócenie czasu trwania symulacji 
bez utraty dokładności [4]. 


4) Warunki brzegowe dla dysku w momencie rozpoczęcia symulacji 


xkt2 = 1. - xkt/1l 
ce = 0. 

fa1(1,1) = 0.0 
fbi(1,1) = 1.0 
li = 0 


gdzie ce jest wartością efektywności zbierania elektrody na początku symu- 
lacji, a 11 jest liczbą wykonanych iteracji. 


5) Blok iteracji 


5A) 


1001 11 = 11 + 1 
test2 = ce 
lxi = amini(float(1x),6. + sqrt(ddmax + 11) ) 


bel określa ilość warstw stanowiących obszar dyfuzji; na początku symulacji 
obszar ten jest mały i rośnie w czasie. Lx1 uwzględnia przyrost warstw w czasie 
aż do osiągnięcia wartości lx odpowiadającej ilości warstw w stanie stacjonarnym. 


5B) DYFUZJA 


do 20 j = 2,1x 
do 20 k = 1,k2 
if(j.eq.1.or.k.eq.1) goto 767 
if(fa1(j,k-1).eq.1..and.fbi(j,k-1).eq.0.0)goto20 
767 fa(j,k) = fai(j,k) + dma + (fai(j+1,k) - 2 * 410. 90 + fai(j-i,k)) 
fb(j,k) = fbi(j,k) + dmb * (fbi(j+1,k) - 2 * fbi(j,k) + fbi(j-1,k)) 
20 continue 


Dwie instrukcje warunkowe pozwalają na uniknięcie niepotrzebnych obli- 
czeń. SC) KONWEKCJA 


do 63 je2, lx1 

inn=in(j) 

do 63 k#1,k2 

1f(j.eq.1.or.k.eq.1) goto 67 

1f(fa1(),k-1).eq.1..and tb t. k-1). eq.0. 0)goto63 

67 irreir(),k) 

if(irr. le. O)goto 64 

af (w()). le. . S. and. wif j,k). le. 5)fai(j,k)efa( inn, irr)+w(J)* 
(e (Inn l. irr)-fa(inn, irr))+w1(j,k)*(fa(inn, irr+i)-falinn, 
eirr)) 

Af (w(j). le. . S. and. wi j,k). 10. 551 (J. k)=fb(inn, irr)+w(j)* 
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(Hot inn- . irr)-fb(inn, irr) tg. K )-ttot inn, irr+1)-fbl(inn, 
*1rr)) 

continue . 

if (w()).gt. .5.and.w1(J,k).gt..5)fa1(),k)efa(inn*1, 1rr+1)-(1-w(J) 
*).(fa(inn*1, irr+1)-fa(inn, trr+1))-(1. w1(j,k))e(fa(inn*l,irr+1)- 
*fa(innel, irr)) 

if (w(3).gt..5.and.w1(j,k).gt. .5)fb1(j,k)efb(innet, 1rr+1)-(1.-w(J) 
*).(fb(inn*1, trr+1)-fbl(inn, 1rr+1))-01.-w1(J,k))e(fblinn*1, irr+1)- 
*fb(inn*1,trr)) 

continue 

if (w(j). le. .5.and.w1(j,k).gt..5)fa1(j,k)=fa(inn, irr+1)+w( j)e(fa( 
*1nn+*1, irr+1)-fa(inn, trr+1))-(1. -wi(),k))*(fa(1nn, irr+1)-fa(inn, 
sirr)) 

if (w(J). le..5.and.wi(J,k).gt. .5)fb1( J,k)=fb(inn, 1rr+1)+w( J)e (tbe 
sinnti, 1rr*1)-fb(1nn, irr*1))-01. w1(),k))e(fbl(inn, irr*1)-fb(lnn, 
eirr)) 

continue 

1 (w(J).gt. .5.and.wi(J,k). le. .S)fal()j,k)=falinnel, trr)=(1. w(j))e( 
*fa(1nn+1,irr)-fa(inn, irr))+wi(j,k)e(fa(inn*l,irr+1)-fa(inn*1,irr)) 

1f(w(J).gt. . S. and. v1 (J. kJ. le. .5)fb1(j,k)=fb(inn+1, irr)-(1. )- 
*fb(1nn+1, irr)-fb(inn, irr))+w1(J,k)*(fb(inn*1,irre1)-fblinn*i,irr)) 

goto63 


64 fa1(J,k)=fa(inn, 1)+w(J)e(fa(inn+1,1)-fa(1nn,1)) 


fb( J, k)=f£b( inn, 1)+w(J)*(fb(inn*1,1)-fbfinn,1)) 


63 continue 


Przy zastępowaniu „starych” stężeń „nowymi” (stężenie w komórce (j,k) jest 
zastępowane stężeniem z komórki (inn,irr)) uwzględniono, że wartości xj i rkk, które 
odzwierciedlają położenie roztworu na początku przedziału At nie pokrywają się ze 
środkiem komórki (inn,irr), w związku z czym wprowadzono poprawki, zakładając 
liniową zmianę stężenia pomiędzy komórkami. Dwie pierwsze instrukcje warunkowe 


pozwalają uniknąć niepotrzebnych obliczeń. 


5C) KINETYKA 


if (xkt.eq.0.0) goto 71 
do 79 j = 1, 1x1 
do 79 k = 1,k2 


70 fb1(j,k) = fbi(j,k) * xkt2 


71 


Jeśli xkt = 0.0, to ten fragment programu nie jest realizowany 


5D) MODYFIKACJA WARUNKÓW BRZEGOWYCH 


fa1(1,1) = 0.0 

fb1(1,1) = fb1(1,1) - dmb + (fbi(1,1) - fb1(2,1)) + dma + fai(2,1) 
do 33 k = 2,ki 

fai(i,k) = fai(i,k) + dma + (fai(2,k) - fai(1,k)) 

fbi(1,k) = fbi(i,k) + dmb = (fbi(2,k) - fbi(1,k)) 

do 32 k = ki+1,k2 

fa1(1,k) = fai(i,k) + dmb * fb1(2,k) - dma + (fai(i,k) - fai(2,k)) 
fbi(1,k) = .0 


5E) OBLICZANIE PRĄDU DYSKU I PRĄDU PIERŚCIENIA 


zri = 0. 
ee = (iri - 6) +2 
do 46 k = ki+1,k2 
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kk = k + ir1i 2 
ra = (((kk + .6) „ 2) - ((kk - .5) + + 2)) % se 
ffb = fbi(2,k) 
45 sri = sri + ffb * ra 
zr = zri + sqrt(dmb) + sqrt(float(1)) 
zdl = rt (daa) + fai(2,1) + sqrt(float(1)) 
ce = zr/zd 


gdzie zr oznacza prąd pierścienia, zdl - prąd dysku, a ce - efektywność zbie- 


5F) SPRAWDZANIE WARUNKÓW KOŃCA SYMULACJI 


rania. 


epsi = fai(2,1) + sqrt(dma*1) + tet 
if (11.gt.20)eps2 = abs(ce-test2) 
if (11.1e.20) goto 1001 

if (abs(eps2-eps).le.0.) goto 1002 
goto 1001 


Parametr epsl oznacza stosunek prądu dysku obliczonego z symulacji 
do prądu dysku w stanie stacjonarnym obliczonego według wzoru Levicha [11]: 
i., = 0.62nFApc?D?2/3w1/2y-1/6, Parametr epsl w idealnym przypadku powinien 
na końcu symulacji wynosić 1. Ponieważ symulacja obarczona jest zawsze jakimś 
błędem odchylenie eps od jedności na końcu procesu wskazuje na dokładność symu- 
lacji. Koniec symulacji następuje gdy różnica pomiędzy efektywnościami zbierania 
w bieżącej i w poprzedniej iteracji jest mniejsza od założonej wartości eps. Należy 
zwrócić uwagę na to aby warunek ten był realizowany dopiero po jakimś czasie 
trwania symulacji, ponieważ na początku procesu efektywność zbierania elektrody 
jest równa lub bliska zeru. Etykieta 1002 odnosi się do części programu związanej 
s operacjami wejścia-wyjścia. Jest to blok kończący symulację. 


WYNIKI 


W tabeli zamieszczonej poniżej zestawiono wartości M (kinetycznej efek- 
tywności sbierania) otrzymane dla elektrody z cienką przerwą izolacyjną i cien- 
kim pierścieniem w wyniku zastosowania przedstawionego programu s wartościami 
otrsymanymi przez Pratera i Barda oras s wartościami otrzymanymi s wzorów ana- 
litycznych 9.8 s [7]. Z porównania tych wartości wynika, śe wartości Nz s przed- 
stawionej symulacji nie odbiegają od wartości otrzymanych przez Pratera i Barda 


Przedstawiony powyżej program został wykorzystany do obliczania stałych 
szybkości reakcji rodanowania pirokarechiny, pirogalolu oraz resorcyny. Reakcje ro- 
danowania powyższych pochodnych fenolu są reakcjami pseudo-pierwszego rzędu. 
Użycie metody symulacji cyfrowej było podyktowane tym, że w badaniach ekspe- 
rymentalnych wykorzystywano elektrody s szerokim pierścieniem i szeroką przerwą 
isolacyjną. Wyniki badań oraz dyskusja zostaną wkrótce opublikowane. 
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Tab. 1. Porównanie wartości kinetycznej efektywności zbierania w procesie ECI 
otrzymanych z symulacji przeprowadzonych wg przedstawionego tu programu, symulacji 
Pratera i Barda oraz ze wzorów analitycznych Albery’ego i in. [7] 


Symulacje wg Wzory 
ME 
programu 
0.1 


analityczne 
0.316 


0.633 
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SUMMARY 


Prater's and Bard's digital simulation technique for investigations EC I process at rotating 
ring-disc electrode has been described in detail. The study presents the manner of treating the 
normal diffusion,normal and radial convection and homogenous kinetics. Moreover, the listing of 
programe in Fortran language and its specification has been attached. 
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